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Maximum Throughput Problem in Dissipative Flow Networks with 
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Abstract — We consider a dissipative flow network that obeys 
the standard linear nodal flow conservation, and where flows 
on edges are driven by potential difference between adjacent 
nodes. We show that in the case when the flow is a monotonically 
increasing function of the potential difference, solution of the 
network flow equations is unique and can be equivalently 
recast as the solution of a strictly convex optimization problem. 
We also analyze the maximum throughput problem on such 
networks seeking to maximize the amount of flow that can be 
delivered to the loads while satisfying bounds on the node poten¬ 
tials. When the dissipation function is differentiable we develop 
a representation of the maximum throughput problem in the 
form of a twice differentiable biconvex optimization problem 
exploiting the variational representation of the network flow 
equations. In the process we prove a special case of a certain 
monotonicity property of dissipative flow networks. When 
the dissipation function follows a power law with exponent 
greater than one, we suggest a mixed Integer convex relaxation 
of the maximum throughput problem. Finally, we illustrate 
application of these general results to balanced, i.e. steady, 
natural gas networks also validating the theory results through 
simulations on a test case. 

I. Introduction 

We consider a network on an undirected graph Q = (V, £), 
where V is the set of nodes with |V| = m and 8 is 
the set of oriented edges with \8\ = n. For a dissipative 
flow network, we have two sets of variables describing the 
network. Each (i,j) S £ is assigned a nominal direction 
i ^ 3 and is characterized by a flow (pij G R from node i 
to node j, and each node i G V is characterized by a node 
potential Wi G R. Such networks which have a description 
via potentials are determined by a set of two equations:“flow 
conservation” and “potential drop relation.” Perhaps the most 
famous example of the dissipative flow networks is a linear 
DC resistive electric circuit, where the flows (pij are currents 
and node potentials are voltage potentials. In this case, 
the flow conservation equations correspond to the so-called 
Kirchoff’s first law, and the potential drop relation is the 
Kirchoff’s second law, or Ohm’s law. Resistive circuits are 
linear dissipative networks, where the potential drop relation 
given by the Ohm’s law is linear. In this paper we consider 
more general dissipative networks where the potential drop 
relations is a general nonlinear relation described by a 
monotonically increasing dissipation function. 

We formulate a variational representation of the network 
flow equations in the form of a strictly convex optimiza- 
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tion problem, thus proving the uniqueness of the solution. 
Additionally, this solution can be obtained efficiently by 
using convex optimization algorithms such as the interior 
point method. The energy function method for gas flow 
networks was studied in [1], [2]. The maximum throughput 
problem is an optimization problem that aims at maximizing 
the amount of flow delivered from sources to loads, while 
satisfying bounds on the potentials and possibly other control 
variables. For general dissipative flow networks, this problem 
is non-convex. In the case when the dissipation function is 
differentiable, we develop a twice differentiable optimization 
formulation of the maximum throughput problem that is 
biconvex in the optimization variables in a certain regime. 
We obtain several properties of the primal and dual form 
of the energy function that are useful to provide gradient 
and hessian information to optimization software during 
numerical implementation. In the process we prove a special 
case of the so-called “monotonicity property” of dissipative 
flow networks, also discussed in a companion paper [3]. 


When the dissipation function takes the form of a power 
law with exponent greater than one, we provide a mixed 
integer convex relaxation of the maximum throughput prob¬ 
lem, where the degree of non-linearity is governed by the 
exponent. Together, we have a framework to obtain a feasible 
solution using the energy function formulation, and then 
check the quality of the solution (degree of sub-optimality) 
using the mixed integer convex formulation. As an example, 
we demonstrate how the setting applies to the natural gas 
networks and also illustrate its performance on a numerical 
example. 


The rest of the manuscript is organized as follows. In 
Section]^ we define dissipative flow networks and equations 
that govern them. Section [^states the maximum throughput 
problem for a dissipative flow network. Section |Vl^ intro¬ 
duces variational formulations of the network flow equations 
and their properties and provides a reformulation of the 
maximum throughput problem using the energy functions. In 
Section we introduce a mixed integer convex relaxation 
for the maximum throughput problem where the dissipation 
function follows a power law with exponent greater than 
one. In Section |Vl] we briefly describe natural gas networks 
and adopt the techniques of the previous Sections to this 
special enabling case. Finally, in Section VII we provide 
computational results on natural gas network models and 


then summarise and present path forward in Section VIII 






II. Dissipative Flow Networks 


A dissipative flow network, where flows on edges are 
driven by potential difference between adjacent nodes, is 
described by the following two sets of the Network Flow 
(NF) equations: 

(a) Flow conservation: 

^ (j)ij = qi, Vz e V (1) 

j&di 

where S R is the injection at node i. 

(b) Potential drop relation: 

- bij, v{z, j) G £. (2) 

The function fij{.) is a differentiable, strictly mono- 
tonically increasing odd function which is called the 
dissipation function. 


Several networks can be described within the framework. 
For example, DC resistive electric circuits [4], AC power 
flows in the lossless approximation [5], [6], [7], [8], balanced 
natural gas networks [9], [10], [2], [11], traffic flows [12], 
[13], [14] etc. The term bij allows us to model additional 
network elements in the actual networks such as compressors 
in gas networks, voltage transformers in DC resistive circuits 
and phase shifters in the approximate AC power flows. 

The system of potential drop relations Q is translation 
invariant, i.e., if {iTi, i G V) is a solution of the NF Eqs. (1 2 1 , 
then a translation of the potentials by a constant, tTi = iti + c, 
is also a solution. To remove this degree of freedom, the 
potential at a particular node (node 1) is fixed at a certain 
value TTi = c. In the case of the resistive circuits, this node 
is called “ground”, in AC power system respective node is 
called “slack bus”, and this term is also adopted to the case 
of the gas networks, where the slack bus (where the pressure 
is fixed) is usually assigned to the largest source of gas in the 
system. Assuming that parameters bij are known, we are now 
left with the number of variables, (| V| + \S\ — 1), equal to the 
number of equations. However, this still does not guarantee 
existence and uniqueness of the NF Eqs. ®. But as we 
will see shortly, in Section II-A this is true for dissipative 
networks with a strictly monotonically increasing dissipation 
function. Further, observe that the set of flow conservation 
equations ([T]) automatically implies total flow balance 


51 7 * = 0 - 


(3) 


iGV 


A. Primal Energy Function and Uniqueness of the Network 
Flow Solution 

Let Fij{x) be the anti-derivative of fij{x), i.e, -^Fij{x) = 
fij{x). The following useful remark follows immediately 
from the fact that fij{x) is monotonically increasing. 

Remark 1: The function Fij{x) is convex. 

The set of NF equations ([T]) and 0 can be rewritten in 
the variational form as the optimal solution of the following 


convex optimization problem [1], [2]: 


min ^ ^ 

(4) 



S.t. ^ ^ (f>ij = (li- 

(5) 

j£di 



The optimization problem in Q-(|^ is clearly convex, since 
by Remark the cost function is convex, and the constraints 
are linear. To see why the solution to Q-Q yields the NF 
equations, we can first form the Lagrangian as 

L{4',Tt)= ^ ^ ^ TTi I ^ ^ fij ~ j > 

(ij)G£ iGV \3&di ) 

(6) 

where tt^ are the dual variables corresponding to the con¬ 
straints 0. To minimize the Lagrangian, we take its deriva¬ 
tive w.r.t. (j) to zero to get 


) — TTi TTj bij ^ (7) 

which means that the dual variables serve the role of the 
node potentials. Since the optimization problem in (Hi-® is 
strictly convex, we also conclude that the NE Eqs. and 
(|^ have a unique solution. 

B. Dual Energy Function Representation 

We now present another variational formulation of the NE 
equations via the dual of the optimization problem in (|^-(|^. 
Since fij{x) is monotonically increasing, it is invertible. Let 
9ijiy) - denote the inverse of fij{x). Let Gij{y) = 

Fifi^j^iy)) - yfr/iy) be the anti-derivative of g^j{y). The 
dual of (|^-(|^ can be expressed in terms of the node potential 
variables as: 

maxy^TT.gi - 5] Gij{TTi - itj + bij). (8) 

iGV (*.i)G£ 

The above problem, when stated as a minimization, is a 
convex optimization problem. Since the Slater’s condition 
[15] always holds for an underdetermined system of linear 
equality constraints, we have strong duality and the optimal 
value of 0-0 and 0 are the same. 

III. The Maximum Throughput Problem 

In a dissipative network the nodes with > 0 are called 
sources and those with < 0 are called sinks. Assume that 
there are constraints on the minimum and maximum value 
of the allowed node potentials, tt^, in the form 

Ti < Ti < Tfi. Vz G V. (9) 

The maximum throughput problem seeks to maximize the 
amount of flow that is delivered to the loads while still 
satisfying the box constraints on node potentials in 0- 
Eormally: 






where 


Maximum Throughput in Dissipative Flow Networks: 


min 5^ CiXi 

(10) 

iGV 

S.t. ^ ^ (pij — Qit Vz G V, 

(11) 

jedi 

TTj TTj — ^ 

(12) 

^ Vi € V. 

(13) 


The potential drop relations ( [T2] ) are in general non-convex 
and make the maximum throughput problem a non-convex 
problem. 


b)].. = -p'(7r, - TTj + h,j), (22) 

j <n-l, {i,j) G £■, 

(23) 

b)] + b^j), (24) 

j&di 

l<i<n-l. (25) 

(c) The inequality 

E{Tr,b) + E*{x,b) — TT^x < 0 (26) 


IV. Representing the maximum throughput 

PROBLEM USING ENERGY FUNCTIONS 

In this Section, we provide another alternative (but, obvi¬ 
ously, equivalent) formulation for the maximum throughput 
problem using the variational representation of the NF equa¬ 
tions in terms of the energy function. To achieve this goal 
we, first, discuss optimal solution and value at the optima of 
the energy function. 

A. Properties of the energy functions 

Define the following convex function given by 

E{n,b)= ^ Gij{'Ki-TTj+bij). (14) 

and let 


has exactly one solution given by tt* = V xE* {x,b). 


Proof: 

(a) The statement in Eq. follows directly from the 
theory of convex conjugates. The proof of ( [20| is 
similar, but we include it here for completeness. Recall 
that E*{x,b) = —E*{b,x). So, it is enough to prove 
that V bF* {b, x) = f*. By definition of (jf, we have 


So, 

dF*{x,b) 

dbki 


4‘*ki + ^ bi 


{iJ)eS 


' bki 




Ml 

bki 


E* {x, b) = sup [tt^x — E[tt, 5)] . 


(15) 


The function E*{x,b) is the convex conjugate of E{TT,b) 
with respect to tt. Analogously, define 


above. 

Remark 2: The function E*{q,b) is equal to the optimal 
value of the primal-dual pair in Q-Q and ([^. Additionally 

F*{b,q) = -E*iq,b). 

Lemma 1: The following statements follow from standard 
results of the convex conjugate and the duality theory [16]. 

(a) Let f* be the unique optimum of Eqs. 0-0 and let 
TT* be the unique optimum of the Eq. 0- Then, 

V,E;*(a:,6) = 7r*, (19) 

VbE*ix,b) = -f*. (20) 

(b) The Hessian of E*{x, b) can be computed with respect 

to X as follows 

VlxE*{x,h)={Vl^E{^,b))-\ (21) 


Using Eq. (j7]l, one derives 

dF*{x,b) 


dbki 


"kl 


bf^i 


(tTj - TTi) , (27) 



F{4,)^ 5] F,M^3)■ 

(16) 



{i,3)eS 



and 





F*{b, q) = sup 51 

(17) 



(Li)e£ 




s.t. 5] = qt. 

(18) 

(b) 


j^di 


(c) 

The 

following remark follows directly from the definitions 
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bki 


(28) 


9<Pti 


Since (l^tj = we have that - 

0. The proof is complete by substituting this result in 
Eq. 


and direct computations. 

By definition —E*{x,b) is the minimum of E{TT,b) 
w.r.t. TT. Since E{Tr, b) is strictly convex, only its unique 
minima tt* can satisfy 


Eor the statement in Lemma[T](b) to be well-defined, we need 
that the matrix V^^i7(7r, b) is invertible. Note that since the 
pressure at node 0 is fixed at ttq = c, indices in Eq. Q 
and Eq. (25 i run from 1 to n — 1. Erom the expressions for 
the elements of the hessian, we derive that V^^i7(7r, b) is a 
(n — 1) X (n— 1) sub-matrix of a weighted graph Laplacian 
matrix with edge weights given by Wij = g'{TTi — tTj + bij). 
Since the graph Q is assumed to be connected, it follows that 
V‘l.„E{TT,b) has full rank and hence is invertible [17], [18]. 

Let us now take a slight detour from the main goal 
(reformulation of the maximum throughput problem) and 
state a certain “monotonicity property” that holds for dis¬ 
sipative EN. Eor a more general statement and proof of the 









monotonicity property without assuming the differentiability 
of the dissipation function fij{<pij) as well as its implication 
on robust optimization on dissipative flow networks, we refer 
to [3], 

Corollary 1: (Monotonicity Property) Consider a dissipa¬ 
tive flow network where the parameters bij and the potential 
TTo are flxed. Let ..., and = 

... Tq^li) be two sets of injections. The injection at 
node 0 is automatically determined by the total flow balance 
Eq. Assume that for all 1 < i < n — 1 one has qi < ri. 
Let = (7rJ^\ .. and = (tt^\ ... be 

the corresponding potentials resulting from the two injection 
profiles. Then one arrives at for all 1 < i < n—1. 


B. Equivalent formulation of the maximum throughput prob¬ 
lem 


Equipped with Lemma [ 

1 one restates the 

maximum 

throughput problem 

(lOi-^ 

13 

' as 


min 

c^x, 



(29) 

s.t. 

E{7T, 

b) 

+ E*{x,b) < 0, 

(30) 


7Li < 

TTi 

< Ti. \/i G V. 

(31) 


This formulation is exactly equivalent to the original 
maximum throughput formulation. Let us present two other 
alternative formulations, which (as we see below) may have 
some practical advantages. 


V. Mixed Integer Convex Relaxation 

In this Section we formulate a mixed integer convex 
relaxation of the maximum throughput problem ([T0|)-([T3]) 
in the special case when the dissipation function f{x) is 
a power law, i.e., f{x) = sgn{x)\x\°- with a > 1. Let us 
start introducing binary variables Sy S {—1,1} associated 
with each edge {i,j) C S that denote the flow direction, i.e., 
Sijfij > 0. We can then restate the non-convex constraint 
( [T^ as 

TTj TTj blj — 1 (34) 

or equivaliently + bif) = ■ (35) 


Next, we apply the McCormick relaxation [19] to the bilinear 
term on the Ihs and relax the equality to an inequality to get 
the following pair of inequality relaxations of ( [T2| ): 

|<^zir < - (TTi - TTj) + {Wi - TTj), (36) 

< Sijin, - Tfj) + (tt, - TTj) - {tt, - Tfj). (37) 


Since a > 1 the constraints ( [36] l and p7| ) are mixed integer 
convex constraints. The McCormick relaxation is tight when 
at least one of the variables is binary, the only relaxation in 
this set of constraints consists in replacing the equality in 
0 by an inequality. 

The mixed integer convex relaxation along with the formu¬ 
lation in Section VI- A| provide the framework for finding a 
feasible solution of the maximum throughput problem as well 
as for testing the quality of the solution through comparison 
with the lower bound. 


VI. Natural Gas Networks 


Energy function formulation 1: 

min c^x, (32) 

s.t. E{TT,b)E*{x,b) < e, 

7Li Ti. \/i e V. 


Energy function formulation 2: 

min (FX-\-M{E{T:,b)-\-E*{x,b) —x) (33) 
s.t. TT^ < TTi < 7fi. Vi € V. 


In (32 1 e > 0 is a small constant. In (33 i M is a large positive 
constant that helps to enforce the constraint ( [30| . 

Eq. @ and Eq. ( [33| are useful because, first, the 
two expressions are twice differentiable and hence allow 
application of second order gradient descent methods for 
local optimization. Second, when b is not an optimization 
variable, the expressions are biconvex in tt and x which 
means that we can use here standard heuristics, such as 
alternating minimization over tt and x. Moreover, solvers 
using second order methods can be provided with first and 
second derivative information by using Eqs. and 

Eqs. (231-(25 1 . The biconvexity is lost when b is also subject 
to optimization, but the twice differentiability remains. 


In this Section, we introduce natural gas networks and 


their properties and apply the techniques of Section VI- 


[A| and |V] to the maximum throughput problem over these 
networks. Natural gas has emerged as an alternative to other 
fuel sources such as coal and fuel oil due to its lower carbon 
footprint. Improved technology for tapping new natural gas 
reserves has increased the availability and reduced its cost 
leading to increasing utilization of natural gas in power 
generation (gas turbines) and household heating. Natural 
gas networks use a system of high pressure transmission 
pipelines to transport gas from sources to loads. The flow of 
gas in pipelines is driven by gradient in pressure along the 
length of the pipe. Pressure falls rapidly with distance along 
the pipelines, so compressor stations are installed to locally 
boost the pressure and maintain the desired pressure and 
throughput at the load end. Typically, compressor stations 
are placed every 50-100 miles. Many of these compressors 
are powered by a fraction of the gas in the transmission 
line itself. The fraction of gas consumed can vary from 
2 — 5%. Due to the increasing demand for natural gas, there 
is an increasing need for optimization of pipeline operational 
planning. In what follows, we describe two such optimization 
problems of interest. 

The first is the problem of minimizing the cost asso¬ 
ciated with compression. In the literature, study of this 
problem can be traced back to [20] (see also [21] for a 













review of this problem) where a Dynamic Programming 
approach was used to solve it. More recently in [22], a 
convex Geometric Programming (GP) based approach was 
developed for compressor cost optimization in tree networks. 
This approach offers several desirable properties such as fast 
convergence, and scope for generalization to a heuristic for 
loopy networks. Since the compressors consume a part of 
the gas itself, the compressor cost optimization problem is 
also related to the next problem below. 

The second is the maximum throughput problem intro¬ 
duced in Section III This is particularly relevant in, for 


example, Norway, which produces more gas than can be do¬ 
mestically consumed. This situation creates strong economic 
incentives to sell the excess gas to the rest of Europe. In 
this case the primary objective of the gas pipeline operator 
becomes to maximize the throughput of the network within 
the engineering and contractual constraints that must be 
satisfied. 

In the rest of this Section, we describe the equations 
governing the flow of natural gas over natural gas networks in 
the balanced (steady state) regime and apply the machinery 
of Section VTA and to this special case. 

Natural gas networks in the steady state are governed by 
the following set of NF equations: 


'y ' 'Pij — 9i) 
j&di 


Pi Pj T ^ij — ^ij4^ij\4'ij\'! 


(38) 

(39) 


where (pij is the flux of gas (flow per unit length) from node 

i to node j in the steady-state, and pi is the pressure at 

node i. The quantity bij models the pressure boost provided 
by the compressor on the edge {i,j) G S. For a more 
detailed derivation of the steady-state gas flow equations 
from dynamic fluid mechanic equations, we refer to [22]. 

From Eqs. P8|39| l it is apparent that natural gas networks 
in steady state can be represented within the framework of 
dissipative flow networks introduced in Section |^by setting 

7ri=Pi, (40) 

and fij{4>i:j) = 5^j(t>^j\(|>ij\■ (41) 


A. Energy function representation for maximum throughput 
in natural gas networks 


In this Section, we adopt the reformulation of the max¬ 
imum throughput problem from Section to the case of the 
natural gas networks. The expressions for functions E and 
F in Eqs. (|14|16[) for gas networks are given by 


£'(7r,6) 


F{f) 


2 

3 

1 

3 


E 


3 

2 


E/ ^^3 I I ^ • 

ii,3)eS 


(42) 

(43) 


The inverse py (y) = is given by 

9ij{y) = sgn{y)^/\^\ (44) 


and the hessian matrix V^^i?(7r,h) is given by 

= - i|7r, - TT^- -f (45) 

{Vl^E{ir,b))^^ = ^3 + ( 46 ) 

j^di 

The energy function formulation of the maximum throughput 
problem for the gas networks can be obtained substituting the 
expressions above to Eqs. ( p9]|3T] i. 

B. Mixed Integer Quadratic Programming (MIQP) relax¬ 
ation of the maximum throughput in the natural gas networks 

The gas network dissipation function ( [39] l follows a power 
law with q: = 2. So the Mixed Integer Convex relaxation 
of Section |V] is now a Mixed Integer Quadratic Program 
(MIQP). We write the full program below for completeness. 


MIQP for Maximum Throughput in natural gas networks 


min 


s.t. 


y^^cjXj 

(47) 

iGV 

M 

II 

H 

(48) 

j^di 

^ SijQti TTj) 

1 

-f 

1 

E 

— Ttj) + 

(TTi - TTj) - (TLi - TTj), 

(50) 

TT < TT < 7f, 

(51) 

X < X < X, 

(52) 

b<b<b. 

(53) 


As remarked before, the only relaxation in the above 
program is the conversion of equality into inequality in 
Eqs. ( 49|50 1 . This corresponds to allowing “decompression”, 
i.e., the pressure drop is now allowed to be larger than 
that implied by the flow. A similar relaxation was used in 
the Geometric Programming approach in [22] to preserve 
convexity, and was observed to perform well in numerical 
experiments. 


VII. Numerical Computations 

In this Section, we test the optimization formulations of 
Sections |VI-A| |V] on a model of natural gas network. Fig. [T] 
illustrates this exemplary network. There are 16 nodes and 
18 edges. The nodes marked with green are sources and 
those with red are sinks. The arrows represent edges and 
the direction of the arrow denotes the nominal orientation of 
the edge. The cost Ci is set to zero for sources and non¬ 
consumers. For consumers, we set q proportional to the 
nominal demand at node i. We consider three different sets 
of upper bounds for the pressure squared variables tt: a) 
Tfi = 10 6) 7fi = 5 c) Tfi = 3.5. In all cases we set tt, = 0.5. 
The slack node 0 is set to the maximum pressure. 

We solve the maximum throughput problem using the 
energy function formulation (|29|)-pT]) and the MIQP formu¬ 
lation ([47]i-([53]l both with and without compressors included 














TABLE I 

Numerical performance for the test without compression 
(See text eor details.) 


Pressure Bounds —> 

0.5-5.0 

0.5-4.25 

0.5-3.5 

Energy function heuristics 

-651 

-576 

-491 

Global Solver 

-665 

-579 

-501 

MIQP 

-685 

-581 

-501 


TABLE II 

Numerical pereormance eor the test with compression (See 

TEXT EOR DETAILS.) 


Pressure Bounds —> 

0.5-5.0 

0.5-4.25 

0.5-3.5 

Energy function heuristics 

-694 

-627 

-584 

Global Solver 

-737 

-632 

-598 

MIQP 

-752 

-640 

-598 



Fig. 1. An exemplary (synthetic) gas network model 


into the set of optimization variables. To test the quality of 
the solution and bounds thus obtained, we compare the result 
with the one obtained with the help of publicly available 
global optimization software. 

We implement the energy function formulation in IPOPT 
[23], which is an interior point method for local optimization. 
We provide IPOPT with first derivative information by using 
Lemma [^a), and allow for automatic numerical Hessian 
computations. For the MIQP implementation, we use BON- 
MIN [24], which is a mixed integer convex programming 
software that uses IPOPT as the local search subroutine. 
We use SCIP [25], a general purpose constraint optimization 
software, for global optimization. 

Table [I] shows the cost comparison between the energy 
function formulation, the MIQP relaxation and SCIP. The 
solution provided by the energy function formulation is 
a feasible solution and hence it is an upper bound. We 
conclude, based on the experiments, that both the lower 
bound obtained by the MIQP and the upper bound from the 
energy function formulation are fairly close to the global 
optimum. 

Table shows the optimal cost achieved by the energy 
function formulation, the MIQP relaxation and by SCIP. 
Similar to the case without compression, one observes that 
the lower and upper bounds found by the MIQP and by the 
energy function based formulation are close to the global 
optimum. 

VIII. Conclusions and Path Forward 

We developed a strictly convex variational representation 
of dissipative network flow equations that proves the unique¬ 
ness of the solution and provides an efficient algorithm to 
find it numerically. For the maximum throughput problem, 
we presented a twice differentiable optimization formulation 
using energy functions as a method to obtain a solution, 
and a mixed integer convex relaxation to test the optimality 
(or degree of sub-optimality) of the solution. We test both 
methods on a model of natural gas networks and observe that 


they were successful in producing solutions close to global 
optimality. 

For future work, it would be useful to extend these 
techniques to dynamic optimization problems on dissipative 
flow networks (such as transient gas networks). The op¬ 
timization formulations could also be used as subroutines 
in multi-level programs associated with network planning 
and expansion. Coupled with the monotonicity property, it 
also provides scope for extension to network optimization 
problems accounting for uncertainty such as the robust 
maximum throughput problem [3] and the problem of finding 
the most probable fluctuations that cause violation of bounds 
known as the instanton problem (see e.g. a similar setting 
discussed in the context of power systems [26]). 
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Appendix 

Proof of Monotonicity Property (Corollary 0-- By 
Lemma [^a), VxE*{x,b) = tt*, which means Va;7r* = 
VxxE*{x,b). By Lemm a [T b). we have the Vl^E*x,b = 
{Wl^E{TT,b)y \ Erom (23 i and (251 we observe that the 
Hessian V‘^^E{TT,b) is a diagonally dominant invertible 
matrix with positive diagonal entries and non-positive off- 
diagonal entries. Erom a standard result in linear algebra we 
have that all entrees in S/^E* {x,b) are non-negative. 
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